﻿/*
 * Krovak
 * 
 * initial implementation : 20.11.2008
 *
 * Copyright (c) Adam Sida, asida@centrum.cz.
 * All rights reserved.
 * 
 * for license see Petra.Gis.License.txt
 * 
 * this is for John Czachurski ;-)
 */
using System;
using System.Collections.Generic;
using System.Runtime.Serialization;
using Petra.Gis.Geometry.Support;
using Petra.Gis.ReferenceSystemServices;


namespace Petra.Gis.ReferenceSystemServices.Transformations
{
  /// <summary>
  /// Implements the Conformal Oblique Conic "double" Krovak projection. 
  /// Bessel to Gauss sphere
  /// Gauss to cone (reversible)
  /// </summary>
  /// <remarks>
  /// <para>The Conformal Oblique Conic Projection of Czechoslovakia was prepared by
  /// survey engineer Josef Krovak, in the year 1922 for construction of cadastral maps
  /// and topographical maps of medium scales for the civil geodetic service of Czechoslovakia.
  /// The projection is based on the Bessel ellipsoid of 1841, which is widely adopted in
  /// Central Europe. Longitude λ is assumed from prime Ferro meridian (an island of the
  /// Canaries), not from Greenwich meridian! The projection is conformal, so that in the
  /// projection plane meridians and parallels intersect at right angles. Angular distance between
  /// Ferro and Greenwich is given as
  /// λ greenwich = λferro - 17º 39' 59.7354"
  /// but round value 17º 40' was used for topographic mapping in Czechia, Slovakia and Austria.
  /// The origin of the S-JTSK coordinate system is located in Finnish Basin, X-axis normally
  /// coincides with the meridian 42º 30' from Ferro with increasing to the South. Y-axis is
  /// perpendicular to the X-axis and increasing to the West.
  /// see more at <see cref="http://krovak.webpark.cz/e_version/krovak.pdf"/>
  /// </para>
  /// </remarks>
  [DataContract(Name="TRANSFORM")]
  public class Krovak : ProjectionTransform
  {

    //private static AxisList ProjectionAxisList;

    private double alfa = 1.000597498372;
    private double u0 = 49.4599572917;
    private double fi0 = 49.5;

    //souradnice posunuteho polu pro definovaneho na gaussove kouli,
    //oproti zemepisnemu polu prenesenemu na kouli
    private double Uk = 59.7118602500 * Math.PI / 180;
    private double Vk = 42.5253936806 * Math.PI / 180;

    private double n = 0.97992470462;
    private double ro0 = 1298039.00462;

    //computation variables
    double kuzpolX, kuzpolY;
    double kouleX, kouleY;
    double kartX, kartY;
    double du, t, u, v, dfy;

    //private double lon_center;		//Center longitude (projection center)
    //private double lat_origin;		//center latitude



    //static Krovak()
    //{

    //}



    /// <summary>
    /// xaml constr.
    /// </summary>
    public Krovak()
      :base()
    {
      //_className = "krovak";
    }




    /// <summary>
    /// Initializes the Krovak object with the specified parameters to project points. 
    /// </summary>
    /// <param name="proj">projected coordinate system with s-jtsk projection</param>
    /// <remarks>
    /// </remarks>
    public Krovak(ProjectedCoordinateSystem proj)
      : this(proj, false) { }






    /// <summary>
    /// Initializes the Krovak object with the specified parameters.
    /// </summary>
    /// <param name="parameters">List of parameters to initialize the projection.</param>
    /// <param name="isInverse">Indicates whether the projection forward (meters to degrees or degrees to meters).</param>
    /// <remarks>
    /// <para>The parameters this projection expects are listed below.</para>
    /// <list type="table">
    /// <listheader><term>Items</term><description>Descriptions</description></listheader>
    /// <item><term>central_meridian</term><description>The longitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the longitude of the point which in the absence of application of false coordinates has grid coordinates of (0,0).</description></item>
    /// <item><term>latitude_of_origin</term><description>The latitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the latitude of the point which in the absence of application of false coordinates has grid coordinates of (0,0).</description></item>
    /// </list>
    /// </remarks>
    public Krovak(ProjectedCoordinateSystem proj, bool isInverse)
      : base(proj, isInverse)
    {
      UpdateTransform();
    }



    ///// <summary>
    ///// s-jtsk axis list
    ///// </summary>
    //public override AxisList Axes
    //{
    //  get
    //  {
    //    // this is default s-jtsk axis list with correct orientation
    //    // this orientation must remain as defined because of projection equations
    //    return new AxisList()
    //    {
    //      new Axis("X", Axis.South),
    //      new Axis("Y", Axis.West)
    //    };
    //  }
    //}







    protected override bool parametersValid()
    {
      //todo:more check
      return true;
    }






    protected override void updateTransform()
    {
      //if (!(Projection.Unit is DistanceUnit))
      //  throw new ArgumentException("krovak projected coordinate system must have DistanceUnit");

      //_projectedUnit = (DistanceUnit)Projection.Unit;

      if (Projection.ProjectionDef.ClassName.ToLower() != "krovak")
        throw new ArgumentException("invalid projected coordinate system, must have Krovak projection");

      Metadata = Projection.ProjectionDef.Metadata;

    }






    /// <summary>
    /// Converts coordinates from lon,lat to projected x,y.
    /// </summary>
    /// <param name="lonlat">The point in radians on bessel ellipsoid</param>
    /// <returns>Point in meters</returns>
    protected override void geographicToProjectedN(ICoordinate lonlat)
    {
      //point kart, koule,kuzpol,Bessel,JTSK
      //Bessel=coord

      //// convert input units to radians
      //lonlat.lon *= _geogUnit.RadiansPerUnit;
      //lonlat.lat *= _geogUnit.RadiansPerUnit;

      //temporary convert to degrees, because all equations below are in degrees ;-(
      lonlat.lon = ProjMath.Radians2Degrees(lonlat.X);
      lonlat.lat = ProjMath.Radians2Degrees(lonlat.Y);

      // Fi, lambda => U, V
      //double alfa = 1.000597498372, u0 = 49.4599572917, fi0 = 49.5;
      dfy = lonlat.lat - fi0;
      kouleX = (u0 + dfy * (99.8585979496 / (100)) + 86.50351075 / (1000000.0) * (dfy * dfy) -
      15.1091 / (100000000.0) * Math.Pow(dfy, 3) - 117.3673 / (10000000000.0) * Math.Pow(dfy, 4)) * Math.PI / 180; //U
      kouleY = Math.PI / 180 * (alfa * lonlat.lon); //V
      // U,V => sirka, delka
      //double Uk= 59.7118602500 * PI / 180;
      //double Vk= 42.5253936806 * PI / 180;
      //Sirka:
      kartX = Math.Asin(Math.Sin(Uk) * Math.Sin(kouleX) + Math.Cos(Uk) * Math.Cos(kouleX) * Math.Cos(Vk - kouleY));
      //Delka:
      kartY = Math.Asin(Math.Sin(Vk - kouleY) * Math.Cos(kouleX) / Math.Cos(kartX));
      // sirku, delku => ro, eps
      //double n = 0.97992470462, ro0 = 1298039.00462;
      //ro:
      kuzpolX = ro0 * (Math.Pow(Math.Tan((39.25 + 45) * Math.PI / 180) / (Math.Tan(kartX / 2 + 45 * Math.PI / 180)), n) );
      kuzpolY = n * kartY; //eps
      // ro, eps => XY
      lonlat.X = kuzpolX * Math.Cos(kuzpolY); //X
      lonlat.Y = kuzpolX * Math.Sin(kuzpolY); //Y

      ////convert to projection units
      //lonlat.X /= _projUnit.MetersPerUnit;
      //lonlat.Y /= _projUnit.MetersPerUnit;

    }






    /// <summary>
    /// Converts coordinates frm projected x,y to lon, lat.
    /// </summary>
    /// <param name="p">Point in meters</param>
    /// <returns>Transformed point in radians on bessel ellipsoid</returns>
    protected override void projectedToGeographicN(ICoordinate coord)
    {
      //sub XY_Bessel(point& P)
      //XY => ro, eps
      //point kart, koule, kuzpol
      //JTSK=P

      //// convert input units to meters
      //coord.X *= _projUnit.MetersPerUnit;
      //coord.Y *= _projUnit.MetersPerUnit;

      kuzpolX = Math.Sqrt(coord.X * coord.X + coord.Y * coord.Y); //ro
      kuzpolY = Math.Atan(coord.Y / coord.X); //eps
      // ro, eps => sirku, delku
      //double n = 0.97992470462, ro0 = 1298039.00462;
      v = Math.Tan((78.5 / 2 + 45) * Math.PI / 180);
      t = ro0 / kuzpolX;
      u = Math.Pow(t, (1 / n));
      kartX = 2*(Math.Atan(v * u)-(45 * Math.PI / 180)); //sirka
      kartY = kuzpolY / n; //delka
      // sirka, delka => U,V
      //double Uk = 59.7118602500 * PI / 180;
      //double Vk = 42.5253936806 * PI / 180;
      kouleX = Math.Asin(Math.Sin(kartX) * Math.Sin(Uk) - Math.Cos(kartX) * Math.Cos(Uk) * Math.Cos(kartY));//U
      kouleY = Vk - Math.Asin(Math.Cos(kartX) * Math.Sin(kartY) / Math.Cos(kouleX)); //V
      // U,V=> fi, lambda
      //double alfa = 1.000597498372, u0=49.4599572917, fi0=49.5;
      du = (kouleX * 180 / Math.PI) - u0;

      // bessel lonlat
      coord.lat = fi0 + 1.001416022789 * du - 86.87150417 / 1000000.0 * (du * du) +
                16.70197 / 100000000.0 * Math.Pow(du,3) + 117.5089 / 10000000000.0 * Math.Pow(du,4); //fi
      coord.lon = ((kouleY * 180 / Math.PI) / alfa); //lambda

      //again temporary convert to degrees, because all equations above are in degrees ;-(
      coord.lon = ProjMath.Degrees2Radians(coord.lon);
      coord.lat = ProjMath.Degrees2Radians(coord.lat);

      //// convert to original geogcs units
      //coord.lon /= _geogUnit.RadiansPerUnit;
      //coord.lat /= _geogUnit.RadiansPerUnit;
    }




    /// <summary>
    /// Returns the inverse of this projection.
    /// </summary>
    /// <returns>IMathTransform that is the reverse of the current projection.</returns>
    public override SRSTransform Inverse()
    {
        return new Krovak(Projection, !IsInverse);
    }





    public override Dimension TransformDimension
    {
      get { return Dimension.dim2D; }
    }





    //public override string[] ProjectionNames
    public static string[] ProjectionNames
    {
      get { return new string[] { "krovak", "s-jtsk", "sjtsk" }; }
    }





    //new public static Krovak Create(ProjectedCoordinateSystem pcs)
    public static SRSTransform Create(ProjectedCoordinateSystem pcs, bool getAsInverse)
    {

      // create empty Krovak
      Krovak kp = new Krovak();

      ProjectedCoordinateSystem tmp;

      // if axes are not same, then replace, but in copy
      // this will keep original coordinate system untouched
      // and avoid error message when create transform
      // with invalid axis definition
      if (!kp.Axes.EqualParams(pcs.Axes, CsCompareType.CompareAllParams))
      {
        tmp = (ProjectedCoordinateSystem)pcs.Clone();
        tmp.Axes = kp.Axes;
      }
      else
        tmp = pcs;

      if (getAsInverse)
        return new Krovak(tmp).Inverse();
      else
        return new Krovak(tmp);

    }


  }
}
